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INTRODUCTION 

In  solTing  an  ordinary  differential  equation,   one  desires  to  find  a 
■aolutlon",  y  ■  F(x),  which  satisfies  the  differential  equation  and  whaterer 
boundary  conditions  are  imposed.     In  ppoblens  whldi  occur  in  practice,  it 
▼ery  frequently  is  not  possible  to  obtain  an  explicit  F(x).    In  some  cases, 
eY«i  nhen  explicit  representations  are  possible,   the  calculation  of  y  for  a 
series  of  ralues  Xj^  invdres  a  prohibitive  amount  of  work  including  usually 
a  series  of  approximations.     Numerical  solutions  involve  the  calculation  of 

values  of  y  for  x  =  Xq,  x x    by  means  of  approximation  methods. 

Usually  the  tabulated  values  of  y  are  desired  for  the  x^'s  which  are,  more 
often  than  not,   equally  spaced.     In  addition,  it  is  customary  to  discuss  the 
aoeoraoy  of  each  derived  approximation. 

Kumrical  mathoda  for  the  solution  of  ordinary  differential  equations 
are  generally  put  into  two  categories:     predictor-corrector  methods  and 
Runge-Kutta  (one  step)  methods.     The  advanUges  of  the  former  methods  are 
their  greater  accuracy  and  error  estimating  ability,  especially  in  systems 
of  any  complexity.     Rnnge-Kutta  methods  have  the  advantage  of  being  self- 
starting  and  easy  to  program  for  the  oosputer.     Neither  of  these  reasons  is 
very  coiq>elling  whan  subroutines  can  be  written  to  handle  systems  of  ordinary 
differ«ntial  equations,  nor  do  they  overcome  their  disadvantages  in  error- 
estimating  ability  and  speed  relative  to  predictor-corrector  methods.     Runge- 
Ktttta  methods,  however,   find  application  in  starting  the  confutation  and  in 
changing  the  interval  size. 

Oonsidering,  then,  Runge-Kutta  methods  for  starting  Uxe  ooiiputation  and, 
perhaps,  dunging  the  interval,  matters  concerning  stability  and  minimisation 
of  ntnd-off  errors  are  not  significant.     Also,  on  modem  ooa^uters,  adnimiea- 


tlon  of  stor«g«  is  beoondng  less  critical.     In  fact,   the  only  criterion  of 
slgnlfioance  in  Judging  Runge-Kutta  methods  in  this  context  Is  minimisation 
of  the  truncation  error.     It  la  the  purpose  of  this  study  to  derive  Runge- 
lutta  methods  of  the  second,   third  and  fourth  orders  whldi  have  minimum 

truncation  error  bounds. 

For  the  derivation,  we  first  consider  the  case  of  integrating  a  single 
first-order  differential  equation  and  later  extend  the  method  to  a  system  of 
n  simultaneous  first-order  equations.  We  then  derive,  for  the  single  first- 
order  differential  equation,  bounds  which  give  the  least  truncation  error. 
It  seetf  x^asonable,  then,  to  assume  that  methods  which  are  best  in  terms  of 
truncation  error  for  one  equation  of  the  first-order  will  be  at  least  nearly 
best  for  most  systems  of  first-order  equations. 

We  are  concerned  with  the  solution  of  the  ^rot-order  dlfferwtial 
•quAtion  problem  given  by 

?    =    f(x,y)  (1) 

dx 

with  starting  value,     yCxg)     «    yg* 

In  addition  we  assume  that  f(x,y)  satisfies  the  conditions  stated  In  the 

following  theorem: 
Theorem  I  -.  If 

i.     f(x,y)  is  defined  and  oontinuotis  in  the  strip 

a  5    X   i  b,        -oo<y<oo        with  a  and  b  finite; 
11.     there  exists  a  constant  L  such  that  for  any  x  <  [a,b] 
and  any  two  numbers  y  and  y*, 


|f(x,y)       -      f(x,y*)j    <       L    |y  .y*| 


i.e.  the  Lipschitz  condition  of  order  1  is  satisfied  for  all  X| 
then,  F(x)  is  continuous  and  differentiable  for  all  x    t    [a,b]    , 
f'(x)  ■  f(x^r(x))  and  F(xq)  =  y^  (i.e.  the  initial  value  problem 
(1)  has  a  unique  solution    7  »  F(x)  for  all  x    c    [a,b]    ). 
We  OHit  the  proof  of  this  theorem.     (See  Henrici    [5.65]   •) 
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aASSICAL  METHODS 
D>ri ration  by  Taylor  Series t 

The  basis  of  all  Runge^utta  nethods  Is  the  expression  of  the  dlfferenoe 
between  the  ralues  of  j  at   Xn«.,   and   x„  by 

t-i 

Where  the  ui^    are  constant,   and     k;  =  >^n  f  ( '<- *  "^^ '*^" ,  y^  ♦  ^ /««j  k^  )   >  (2*) 

with    «H,»  0  ,  ^,'  Xn..  -  x„    ,     Given  the  weights  w;  ,    the  parameters     «Kt  ,  a^. 
and  using  equation  (2)  we  can  scire  equation  (1). 

We  desire  to  determine  the  w^ ,   m;    ,  a.    so  that  the  coefficients  of 
in  the  Taylor  series  expansion  of  both  sides  of  equation  (2)  about  the  point 
(Xjj.jjj)  are  identical  for  r  a  1,2...,m  for  soner  fixed  m. 

ftqpandlng  F(x^^^)  we  hare, 

2.1 
Hence,  Y**'"y*   *    r(x«».)  -  Rx^) 

S      (Xn».-Xn)rU«)        -^      (X>«.-X»r    r    (Xn\  ♦     ...     . 

iT" 
Bat,       •  m^  •    x,»,  -  Xn   , 

•o  that 

2'.  3.' 


or  aore  generally, 

>••    y-    '  irr-      tr     •  (3) 


Ib  addlUon,  sine*      f(x.»,y«)=  y^       ,  ftssuiilng  diff«r«tlablllty,  it 
fbUova  that 


and 


(*) 


nhara 


^  -  ((«.^)  sf^y 


Thus  wa  haTa  , 


yn»«  ~  Yt« 


*•. 


^(lir*f}yr^<M"> 


(5) 


Mow  daflna  D  ■  ll  ♦  (.  i_ 


where   f,=  ((x«.y-)  ; 


than 


(f.*fiyn(M-)l  -  ^^f*/,w 


(6) 


(7) 


Zb  gaaaral. 


D"  = 


fc<* 


D(D"u)  =  D"*:.  ^^mo^'iy   , 


D'm.  «  w.   . 


Th»  •xpmslon  of  equation  (5)  gl"v»s» 

r--y- -^f- 4^,(1. -f^)''  ^(^-f^n  *■■■' 

nd,  using  •qUAtlon  (7),  va  have 

5  (8) 

♦  7f,Dfog][  -  oat) . 

To  obtain  the  expansion  of  the  right  hand  side  of  equation  (2),  we  use 
the  Taylor  series  expansion  for  tvo  rarlables    [l6,227] 


or 


Let  D^  ■  (-1  j;  ♦(Z:^.jU,i_)      J  on  factoring  out  the  X  •      ^^0) 
ve  have 

U...-.i...y..(2:^jU.f.l  =  Zl^^^|^  •       (11) 

Hev,  using  equations  (10) and  (11),  ve  can  expand  each  ki  on  the  right 
hand  side  of  equation  (2)) 


(12) 


(13) 


|<,=   J«,f(x,  +  -,J%„,    Y"*  ?^/«'jM  > 
or,  g«n«raliting,  m  havv. 


k.=  J^.|ll^-'''-|^^-/V^W%^    ■ 


(15) 


Q«lng  th«  rvaoLta  for  k  ,  j « «^     ,  to  ml  to  ki  «8  an  «qpression  in  povors  of 
iin  ,  wo  haT», 


(16) 


ki  z     Jijcx^.y.)    *■  V„  DJ(/..y.)     ^    ix,  DCfCXn.V,) 


i! 


k,' 


Jj  4'.  '     ' 

2.!  3!  /^. 

aiallarly,  tho  dorlred  equations  for  ko  and  k^,  are. 


(17) 


(18) 


-  „,-  .■..-j,..fl;»-5„,;,.,.^.^^.j^-s.|.-»-)p(  .1 


■T 
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•  i/jnMD;i,  ^-kfl.AJDZfy  -Ay^Jyi^jDjj    ♦  oa\)  . 

Equations  (16)  through  (19)  will  enable  us  to  derslop  all  Runge-Kutta 
Mthods  through  order  four,  and  the  terms  Involving  A^  will  facilitate  the 
diseossion  of  the  error  texms. 

SviMtltating  the  expressions  from  equations  (I6)  through  (19)  and  equa- 
tion (6)  into  equation  (2)  and  equating  powers  of  K  through  X^    we  have, 

U>.  ♦  uJv  ♦u),*«*>4.=.|,  ^20) 

.-.   u,,DJ  *  w^D,;  4  w^O^I  =  Ti  '  (21) 

(22) 


=  ilD*f*f,D^^  *3DfD(^  w;dij 


(23) 


aiaoe  the  <Xi  ,  ^^  and  w;  are  to  be  independent  of  f(x,7),  equations  (20) 
thnagh  (23)  aotually  represent  eight  equations.  Therefore,  the  expressions 


x  " 


in  braokvta  In  th«s«  •quatlons,  which  are  homogeneous  In  the  operators,  must 
equal  the  corresponding  terms  on  the  right.  ia.so  If  these  eight  equations 
are  to  be  Independent  of  f(x,y)  then  the  ratios, 

■ost  be  constant.  This  vUl  be  true  If 

'  ,     '  ■  "  -*, 

1- « 

•««.    =    X         &n  ,  i«    1,3,4     .  (25) 

Tterefore, 

D;   ^    •'iD    .  (26)     , 

finally,  these  eight  equations  become. 


3 


«»«*x/»jt«-     w*('<i/i*i   ♦  -t^^**)         =     -^      > 


U»,#t»    ♦     U|j.<,        *       ^K  -t^  =       -3-       ♦ 


4 


(27) 


In  the  equations  abore,  the  first  corresponds  to  equation  (20),  the  second 
te  equation  (21),  the  third  and  foiurth  to  equation  (22)  and  the  last  four  to 
equation  (23). 

Cooslderlng  then,  equations  (25)  and  (27),  we  have  eleven  equations  In 
13  unknowns,  which  wHI  generallj  be  sufficient  to  determine  the  parameters 
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on  assigning  Tsluas  to  the  two  free  parameters*    ¥s  will  now  consider  the 
Ronge^utta  Mthods  of  orders  2,   3,   and  k, 

Fbr  the  case  when  ■  «  2,  the  system  of  equations  (2?)  retains  only  the 
equations  pertaining  to   J^\   .    These,  with  equation  (25)  for  1  «  2  are: 

<x^^   '    Y    -  (28) 

Three  second  order  equations  of  Interest  correspond  to  the  following  raluesi 

o<2     ,     1/2,  2/3.    1     . 
Substituting  in  equation  (1)  we  have,  respectively, 

)n^-y*    «  i.Hxn*^^'",  ^.^kX>J^)   ,  (29) 


(30) 
(31) 


7-». -yn    »    -^J^nt^Cx^.y.)    +    ^(<«*i.  ,  y„*  in^O]    . 
It  is  Interesting  to  note  in  the  above  equations  that  equation  (29)  is 
the  Vewton.Cotes  open  type  formula  (or  the  improved  polygon  method  or  the 
modified  Euler  mothod)  and  if  f(x,y)  is  a  function  of  x  only,  it  reduces  to 
the  mid-point  rule  of  numerical  integration.     Also,   equation  (31)  is  like  the 
familiar  trapesoidal  rule  when  f(x,y)  is  a  function  of  x  only,  otherwise  it 
is  the  familiar  Runge.Kutta  seoond^rder  method  (or  the  Improved  Suler  method 
•r  the  Hem  method)  most  oomooonly  seen. 


Tot  th«  oasf  «h«n  ■  «  3t  v«  have. 


U),    ♦      CO  ^    4     vO  ^ 


11 


•<j»»)^    ♦     o<»  «0^    =     -      t 


<*W^-»      .<*to,        =    -^ 


(32) 


o<  J  -       fin    •*■  fitt. 

This  is  A  two  p«r«Mt«r  fanily  which  can  be  rewritten  as. 


to, 


'I  = 


3m,  -  Z 


W^  r 


Z-   3««, 


Um»(.<,-^») 


(33) 


a  <, 


fi.     » 


/-M  - 


where   •«%  "^  •«,   ,•<«.•<»**,  •<*  *  \ 


-ia-j_f^»-  •■« "Tt""^":""* 
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»d. 


Two  ooHBon  third  order  83r8teiB8  arst 

kt  =  >t.J(x„*^^. ,  7«*ik.)  , 


(34) 


,.  %: 


(35) 


Iqtution  (35)  aboTB  is  siallar  to  tho  ooaoaon  Simpson's  nil*  i^en  f(x,7)  Is  « 
fteetion  of  X  onlj. 


■  .1;.:.  /■■ 
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For  th«  o«s«  when  ■  «  'f,  we  again  have  a  two  paraiMitar  faoUy  of  oqvu- 


tioos  vhioh  can  be  solred  to  give  the  following  t 


VU,  a     A_   ♦     


tJ,    >  •^**»   -  ' 


^      \iM»c*,..<,.X'--*») 


(36) 


2iK,  (I-  ^oil.) 


/**   '       l«ti(«*»-«*x)[t.«,.,,  -  41H».».<,)  4.3] 

where    «Xt,-<,  #o  ,  o^t  ,»<    :^  I  ^    •<i^-<,      and  the 
denominators  of  the    yd's  do  not  vanish. 
Ihe  Bost  ooHK>n  fourth  order  method  is  obtained  when  we  let    '^  *    '^  *  1/2, 

where         U,  -    >^.sV(yn,^«)     , 

kj  -  >JCx**^^Jk«,  Y.*iU,)    , 


■  ), 


> 


1<» 


Th«  U8«  of  P>de  Approximanta  In  the  Derivation  of  the  Second-Ordar  Method: 

We  will  dlseuss  the  use  of  Fade  approxlnants  In  the  oonstructlon  of 
difference  equations  which  wUl  lead  to  a  fanlllar  numerical  solution  of 
differential  equations*     As  before,  we  are  giren  the  following 

^    a    f(x,y).    with    y(xQ)    =    Jq.  (39) 

We  define  Fade  i^roxinants  4fter  the  following  discussion  which  noti- 
rates  the  definition:  Suppose  we  hare  a  function  P(x)  for  irtdoh  the  Taylor 
series  can  be  written, 

P(X)  =    c.  ♦    t.x   ♦    c^x»   +    ..  .        =    Z—   C^K**      . 

k  »  • 

Ut 


Then  D  (x)  is  a  polynonlal  of  degree  p  in  x.  We  font  the  product  r(x)  L^(k) 
whidk,  after  siifjllfl cation,  yields: 

We  hare  (p  -f  1)  parameters  t.  ,  (k  =  0, 1,  .  .  ,  ,p)  which  can  be  chosen 
so  that  the  coefficients  of  X***""  will  vanish  identically,  for  r  »  1,2,  ,,, 
,p.  Let  V  (x)  be  the  polynomial  of  degree  less  than  (q  -f  1)  fornsd  hy  ell 
terms  of  degree  less  than  (q  <«■  1)  on  the  right  hand  side  of  equatim  ClO). 
Then, 


••/ 
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Th«  rational  fraction,  f(x)  =  i^   '       Is  tho  [0,0,]  Pad©  approxlnant  to  th« 
Anotlon  P(x),  vhere  In  the  notation  [p,a-J  *  P  denotes  the  degree  of  the 
denosLnator  and  q  denotes  the  degree  of  the  numerator. 

This  Fade  approxlmant  enjoys  the  following  basic  properties: 

1.  It  Is  a  uniquely  detemlned  rational  fraction  approximation 
to  P(x). 

2.  If  the  Fade  approxlnant  Cp,^]  is  expanded  In  a  Taylor  series, 
the  first  (p»q-*-l)  terns  will  be  Identical  with  the  first  (p-tq-fl)  terms  In  the 
Taylor  series  expansion  of  the  original  function. 

3*  An  estimate  of  the  error  InvolTsd  In  a  glTsn  Fade  approxlmant  . 
can  be  caLculated  from  the  remainder  term  In  the  Taylor  series  expansion. 

Kopal  [9.163]  regards  property  2.  as  being  fundamental  and  from  It  is 
seen  that  the  first  (p4q-f1)  terms  of  the  Taylor  series  expansion  form  a  Fade 
approxlmant,  namely  the  [o,  p4ql  approxlmant. 

It  can  be  shown  that  the  Fade  approxlmants  [p.p]   to  ln(l  -  x)  eralu- 

ated  near  x  *  -1  are  more  acciirate  than  the  [o,  2p1  Taylor  series. 

We  denote  z.    by  the  operator  D,  and  h  as  the  IntenraL  step  size;  then 
dx 

3^  ■  Xq  -f  hi  ,  yj^  «  7{x^)»    Define  the  forward  and  backward  differences, 
respeetlTily  as  follows: 

A  y^  -  y^^^  -  ^1  • 
V  yi  «  ^  -  y^.l  . 
Using  this  notation,  we  can  derlTa  the  following  symbolic  relationship: 


(*i) 
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Now,  using  the  operator,  as  defined.  In  the  original  equation  (39)  end 
•ubttitutlng  this  result  in  equation  Cfl),  we  have, 

^    «    I)3r    «    f(x,y)     , 

dx  '    .  • 


ilr^jy  ^  f<M)  •       :  («) 


For  a  partioolar  x  and  7,  equation  ClS)  beoones 

where         h  =  t<^iy*) 

Z(yi  -y,..)   =    Zl\fi  -   i»  ^f;      , 

Ib  teras  of  t^,  we  hare 

7.  =  y...  .J,  [f.-i  ((,-(...)]  , 


or 


(^3) 


Iquation  (^3)  la  precisely  the  siBfxLified  Runge^utta  seoond-order  femnla 
with  m  •rror   OW, 
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0«o— trie  IntarpretAtlont 

Oonsldtr  again  th«  dlfferontlal  aquation, 
^  =  ■f(x.y)       ,     y(X.)  -.  y.     , 

and  usa  tha  ordinary  fourth  order  nothod. 


(^) 


w   » 


vharo. 


Tha  gaonatrio  slgnlfloanea  of  thesa  formulae  is  sboim  In  Figure  1, 
There  it  is  aeaa  that 


(*5) 


k^.y) 


fb 


la  the  slope  at    R    of  the  eunre  y  »  F(x)  vhldi  satisfies  the  differential 

equation  in  equation  (<»4).  Thus  from  equation('f5)we  have, 

J|^    i    H*-.>f-)  =   the  slope  of  the  curve  y  s  f(x)  at  (3^tyn)f 

X  "  *•    ''  =    the  slope  of  the  approximating  function  at 

—.  s  H*«'»-r>«,  V**  tk»y  s   the  slope  of  the  approximating  function  at 

k   p 

-^  t   Ric.»i^  «Y««  k,)  s  the  slope  of  the  approximating  function  at 


IB 


Figure  1 


lach  of  thM*  straight  lines  whan  drawn  at  (^fTj.)  Ulustratas  ths  signifi. 
oaaoa  of  ths  corresponding  k^{  1  «  1,  2,  3,   4. 
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TRUNCATION  ERROR 

As  WAS  st«t«d  b«for«,  when  using  the  Runge-Eutta  method  for  starting  a 

solntlon  and/or  changing  the  Interval,  matters  such  as  stability  (propagated 

•rroT)t  and  roundoff  error,  are  not  significant.     Vfo  would  like  then  to 

■Lnlxlae,   as  best  we  can,   the  truncation  (or  discretisation)  error.     To  do 

this  we  first  look  at  the  error  terms  for  the  three  methods  msntlcmed  earlier. 

Kquatlon  (2)  Is  to  be  exact  for  powers  of  h_,  19  to  h*;  hence  the  trunoa- 

n        n 

tion  error  T  oan  be  expressed  as, 

X  =  Y.Xr  *  OUT)   ,  (1,6)  , 

vhora  both  Ym  and  T^  are  dependent  vq>on  the  function  f(x,7)*  In  order  to 
estimate  Tn,  we  hare  to  consider  only  Y  ^^  since  the  bounds  that  we  will  glTS 
K  ^  are  so  conservative  (i.e.  the  true  magnitude  of  Y^j  will  generally  be 
much  less  than  the  bound  we  give)  that  the  term  OCi^n  /will  be  vejy  small 
compared  with  Y^hJ*^  (which  should  be  the  case  if  h^  is  small).  Then  the 
bound  on  Vj^   will  tisually  boiind  the  entire  error  term. 

Using  equations  (8),  (16)  to  (19)  and  (26)  we  can  calculate  the  terms 
^2t  Vj  «»d  Yj^'  'or  m  «  2,  from  equations  (16)  to  (19)  and  equation  (8) 
we  have,  ,  ^ 

n.  =  i^(D>{ .  i,Df)-  w.JL,o:f , 

*»  Z  ' 

and,  usinc  •quation  (26), 


Henoe, 


if?) 
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aiKlI«r  rssults  o«n  b«  derl-ved  for     V3  and    Xk'    ^"  algebraic  manipula- 
tion, howarar,   is  too  tedious  to  ba  given  here.     The  reeuLts  ares 


•V  =     V  •«-«> 


U3v  M*   *     i^x<*     -»     U>4e<4' 


i4 


•)d*^ 


(rf. 


-     ^»»<t»< 


U»*«<4 


-       »*> 


(rb- 

^  >  to 
(JL.   -    ^ 


(^) 


In  order  now  to  bound   V^,  we  assume  the  following  bounds  on  f(x«j) 
its  deriratiTM  in  a  nelghboiiiood  of  (3'^.y„)t 

|f(x,y)|  <   M  , 


(50) 


where  i  4  J   ^    M  and  M  and  L  are  constants,  such  that  M  2  1. 


21 


H«iM  th«  bounds  for  Y,  can  be  deterodned  by, 


\-.  (-t-  ^')dV  '  t^,W    . 


Kov  froa  oquation  (6), 


»»•  >»>/  if 


and,  vBlag  •quatlon  (50), 

> 

^.     <     i^      , 

11^  <^ 

» 

|W|   ^      2.ML 

SUUarly, 

» 

rl 
>^- 

M 

• 

Hanea, 


|D*f|    ^    H12    ♦    2.  ML"    *   ML'     =    ^-ML*-    , 


and 


IY.I    ^(4|i-    •'■^'l    .-i-)fAL^ 


(51) 


In  a  aiallar  vaj,  wa  can  deteradne  bounds  for  Yo  and  Yu,    Tha  derira- 
tloB,  bovarar,  will  ba  osLttad* 


r,V  -r.*-'  *•» 


"•«".'.''  r  - ".' 


:;^  ";:-^j?r^'T?i^^;S33?'^T^5?i5:i 
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(52) 


(53) 


I    /   .» 


^.  =  ]i;  -  •^-<'-'^3)/».Y. 


*s  ♦"♦  , 


* 
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R»f»rrlng  to  equAtions  (2?),   (28)  and  (33),   for  which  the  paraneters  ar« 
vndardetendned,  In  all  three  cases  the  underdetemdned  parameters  will  be 
assigned  ralues  so  as  to  ndninize  the  bounds  we  hars  Jtupt  set  on  the    ]C^» 

Mow  let  us  consider  the  second  order  system  (i.e.  ■  »  2)  and  the  one 
paraMter  family  in  equation  (28)  resulting  from  equating  like  powers  of  h,^ 
up    through  \  * 


We  h*T»  thm,         ^ 


=  /J. 


to.  -     \  —  — *       » 


icK 


X 


W.-     i;^^  .  V      (5») 


»t.    -     -tx 


Siaoe, 


«•  hjiTe,  after  substituting      ^*-'  J^       » 


i55) 


It  is  dear  from  relation  (55)  that  if    ^2  *  2/3,      Y,  is  minimized  and, 


reoTer,    Vg  is  strioUy  less  than  -t"  .     It  is  seen  that  for     *^2  "  2/3  *• 


'2 
hAve  the  following, 

vfaieh  was  one  of  the  second  order  systems  giran  in  equation  (30)  • 


Zk 


Ibr  ■  «  3.  oonslder  the  two  parameter  family  giTwi  In  equatlcKi  (33)  • 
After  solTing  for  w,   and  w^  in  tenw  of  Tf<j,  let    «»<  ^  «  0,     o<  ^  «  2/3  and 
obtain  the  following  one  paraneter  family  of  equations. 


'^x  - 

X 

"3       ' 

h-^ 

U),      = 

±-    -U.,    , 

(56) 


U)1 


Siailarly  we  oan  let  •^g  "  '^3*2/3  with  the  resulting  one  paraaeter  fasily 
ef  equations  being  giraa  by 


=^  1-  -  oJx  ,  (57) 


U>^  *   ^  _  W,   , 


/*»^  ^  4u 


» 

5 


We  again  follow  the  mthod  used  for  n  s  2  and  see  that  the  coefficient 
of  YL^  in  relation  (52)  will  be  a  Minimum  if  0^2  "  ^/^»  ^^3  "  3/^,  which 

giTSS, 

(58) 
Using  the  ralues  girm  in  the  system  of  eqxxations  (56)  we  hsTB  the  following 


bound  for     Y^* 


U3I  ^3^^'     '  (59) 


or,  by  using  those  values  in  the  system  of  equations  {5?)  "^f^  hsTS  as  a  bound 


-  'r 


U.\   <    iHL^   . 


(60) 
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Thus  Illation  (58)  obviously  will  give  the  least  bound  on  the  truncation 
error*  Using  the  raLues  of  the  free  paraneters  (i.e.  '><'^   «  l/2,  *^3  "  3/**) 
nhidi  vere  used  to  arrive  at  relation  (58)  w  have  the  following  third  order 
equaUon  which  will  yield  vidues  of  y  which  are  best  In  the  'least*  trunca- 
tion error  sense, 

where    U,  «  vlvn^Cx^.y^")  ,  (61) 

k,  =  i4(>^"*  \M.  -.  V-*  ^  W)  . 
For  ■  ■  4  we  have  the  two  parameter  family  given  in  equation  (37)  irtiose 

possible  solutions  are  listed  below  for  various  (common)  dioioes  of  the  free 

paraaiters. 

*  '      ■'  '  4      .   ■ 

u,^-     ^  -  ^,  ,  ^,,  ^    ^      ,    .  (62) 


I  ' 


•^1  »    •<^  '     »  ,  -Cs  ^    -^        , 


*^'  -  ■<:     >    u}x  -  ^  -  ^^ 


«-»*  ^     >      /»3t  --  "k 


yfl**    -    -  t      fi^^    ~     — — 


(63) 
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Using  th«  bound  on     Y  |.  as  given  in  relation  (53)  i»ith  some  lengthy 
oo^MUtion,  it  can  be  shovn  that  the  coeffLaient  of  Ifi'^  will  be  ainiBixed 
vhm     ^2  -  0.4  ,     «»<^  a  7/8  -  3/l6Tr  . 

We  hare  then, 

lY^l    <     5.4-(.   >*  10"^  IM""    .  (65) 


Ve  ean  ooapare  this  bound  with  others  for  various  valtiss  of  the  coefficients. 
In  the  system  of  equations  (62)  with  w^  «  5/8,  we  hare, 

U^l  <  1.11  *   10-^  hU^  .  .:^  };\   (66) 

In  the  systMi  of  equations  (63)  with  Wj|^  =  IO/51,  we  hare, 

|)(^)  <  \S.T1  X  lO""-  ML*-  .  (67) 

la  the  system  of  equations  (64)  with  w^  s  -5/78,  we  have, 

lUl  <  n.(.A-  '^iQ-^hL'-  . 

A  mdx  earlier  bound  was  found  by  Lotkin  [9.130]  for  the  system  given  in 
the  system  of  equations  (62)  with  w  s  1/3  and  was 

UA  <    10.14-  '^  »o-^rAL^    . 

fbr  the  sake  of  illustration,  we  will  list  some  of  the  more  oomnon  error 
bovids  previously  derived  for  the  fourth  order  method.    The  bound  for  the 
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dastioal  fourth  order  msthod  with  <^2  "  ^'^^  *^3  =*  O*^  ^s  glTen  by, 
Fbr  tho  Bithod  by  KutU  with  «^2  "  ^Z^,  ^K^  =  2/3  wo  haro, 

kxd   fbr  tho  BBthod  by  0111.  with  ^Kg  «  l/2,  «^^  «  1/2,  w^  «  1  +  -  wo  h«TO, 

U^l  <.  8.83  >*  10-"  ML*-  ,      r^'^®^]  • 

Itt  all  cases  it  is  obvious  that  the  BdniniiM  error  bound  is  giTon  by 
relation  (65).     This  was  first  derived  by  Ralston     [lO,435l    •     Using  the 
ralues  of  the  free  parasneters  which  were  used  to  arrive  at  relation  (65)  we 
have  the  following  minimua  fourth  order  Runge^utta  Bwthodt 

y^^  -  y    -  .17'>76028  k^  -  .551'*8066  kg  +  1.20553560  k^  +  .17118^8  k,^, 

where 

"    ■       "  -■  .   .,r 

:■,-■* 

,■■"'■  ,  ■  '' 

kj  =    ii  J  (x„  +  .4-55  737  i5  i\„  ,  (68) 

y„*   .29(»«177(.I  k,    ♦    .  IS875«^  (,4-k^)  ^  /' 

* 

Bere  y  . .  is  n  approximation  to  the  solution  of  the  given  differential 
n+i 

•qvatioa*  ■'::j'~-:""':.y  ".^  -■ 

I-  :  '  •   •      ,  ■   •  ■         .    '.   ■■  '--  '        ' 

'    -  * 

■  ■    ■  .'■■-■.  ■■-.■:'■''  ,' '  "•  ■..  . 
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THE  RUNQEJttJTTA  METHOD  AS  APPLIED  TO  A  SYSTEM  OF 
FIRST  ORDER  DrFFERDJTIAL  EQUATIONS 

X  87Si«m  of  oirlLnary  differential  equations  of  the  first  ordor  Is  a 
•ystoa  of  oquatlons  of  the  form 

f  =  {'(".V'.y' Y')  , 

r-n^r-r )">>  m 

•     •     a 

irtier*  f  ,f  \  ...,^       *r«  given  functions  of  (a+1)  argumants,  and  the  ssoond 
supersczipt  on  the  y's  refers  to  the  derivative  with  respect  to  x. 

A  set  of  functions    >^H>«),^V*)  ,  ...  ,  y'^J^)      whi<di  are  defined  and 
dlfferentLahLe  in  an  interval  [a.b]    and  satisfy  identioaliy  in  x  the  ^ 

relation,  -. 

Y^'  -  t'<x,>/'U)  .yv*),  ....  Y  *''<))  ,  (70) 

is  called  a  solution  to  the  system. 

The  problem  whidi  arises  frequently  in  practice  and  which  we  will  dis- 
eoss  here  is  to  find  a  solution  of  the  system  of  equations  (69)  irtilch  satis- 
fies the  iBlUal  conditions 

7*^0  =   >);       ,  i=  «,u,  ...  ,  s     ,  (71) 

where  the    *h   are  prsassigned  constants* 

It  should  be  noted  In  passing  that  other  conditions,  more  complicated 
than  these  of  equation  (71),  else  arise  in  practice,  but  will  not  be  discussed 
here* 
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SysUas  of  ordinary  differontlal  aquations  arise  in  several  ways;  two 
general  situations  are  given  in  the  following  pages: 

(i)  Theoretically,  every  ordinary  differential  equation  of  order  higher 
than  the  first  can  be  reduced  to  a  system  of  first  order  equations,  Gbnsider 
the  differential  equation  of  order  a  given  by 

T'-U^.i.r.r-  •■  .7'-"),  ■     *^' 

whei^  f  is  a  given  function  of  (■*!)  argUBonts.  The  reduction  of  this  equa. 
tioa  to  a  system  of  equations  (69)  is  aooonpLished  by  setting 

lev  if  the  fttnetions  V^'7^''"  ^V"*  satisfy  the  systen  of  equations, 

^  (73) 

•  •  • 

7""  '  U^.t.y\  ...  ,Y")  , 

then  the  fmction  y(x)  «  y^x)  will  satisfy  equation  (72).    Thus  the  system 
of  equations  (73)  is  a  special  case  of  the  system  of  equations  (69)  • 

Some  authorities  (MLLne    [l2,82]     and  GiU^«9$    reooammid  this  reduc 
tien  of  hif^er  order  equations  to  a  system  of  equations  of  the  first  order 
also  for  nxaerical  purposes;  others  (such  as  Collatz   [2, 1171     )take  the  oppo- 
site position,   arguing  that  reduction  to  a  first  order  system  increases  boUi 
the  en\>r  nd  the  necessaxy  number  of  operations.     Methods  for  the  direct 
integration  of  equations  of  higher  order  can  be  found  in  most  advanced  texts 
en  the  subject.     However,  the  theoretical  and  experimental  results  presented 
dealing  with  one  equation  will  be  at  least  nearly  best  for  most  systems  of     ■''' 
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•qtutlons.     Bvidanc*  indicates  th*t  it  is  possible  to  control  the  truncation 
error  and  it  can  be  shovn  that  the  rovmd-off  error  is  frequently  substantially 
decreased  when  an  equation  of  higher  order  is  first  reduced  to  a  first  order 
•ystea  and  then  solred  by  an  equivalent  method  for  such  a  system.     See 

Henrici     [6. 123]   . 

(ii)    Systems  of  ordinary  differential  equations  also  arise  in  a  natural 
Way  from  many  physical  problems.     Classical  examples  are  electric  circuits 
with  moz^  than  one  loop  and  mechanical  problems  with  several  degrees  of  free- 
dom.    More  speoifio  examples  are  the  equations  of  motion  of  a  gyroscope,  the 
ftndamental  equations  of  exterior  ballistics  and  the  equations  governing  the 
flights  of  rockets  and  missiles. 

Vector  Notation : 


It  will  be  convenient  for  us  at  this  time  to  simplify  the  subsequent 
analysis  both  oonoeptually  and  formally  by  considering  the  quantities 
y^,  i  ■  1,2, •••   ,s     as  components  of  the  vector 


ConseqwnUy  we  write      ^\%,>j\)\...  .>/*)     ~  f\x,^) 
ComblBe  the  s  fonotions  of    Tc^^.y)  into  another  vector: 


fc'^.y) 


\\\kV 


t  :.. 


If,  \'; 
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Writ*  •quatlon  (69)  in  the  nav  ooiqMOt  fom  as 


Y*  -Tcii.Y) 


(7k) 


If  «•  d«fiiM  th«  Tvotor 


1 


by 


V 
tb«i  ih«  iBiUal  oondlUon  In  (3)  is  gi^^wi  bgr. 


(75) 


la  addition,  as  in  tha  casa  of  one  aquation,  wa  assuM  tha  following, 
that  tha  raotor.Talued  f\mctlon  f(x,y)  of  the  scalar  rariabla  x  and  tha 
Taotor  Y*  (Y',Y*,-..y*)  satisfy  tha  following  two  hypotheses  J 

(i)  ^(x,^)   l8  defined  and  continuous  in  tha  region 


^ii'=^;r 


a   £    X  £   b,    -oo<y*-< 


oo 


i  «  1,2. 


•    •   •   fi 


(ii)     there  exists  a  constant  L  such  that  for  sona  x  «  [a,b]       and  any 
two  reetors   ^  "»•*  Y*  • 

•    l|?<Mi  -?(v-r)ll  ^  Uh-vi  . 

where  Hvll   iadleatas' the  norm  of  the  rector    v  . 

Wt  ndt  the  lengthy  proof  of  the  following  theoren     \/f,113l  • 
Theoreat    Let  the  function    ?^X,^)      satisfy  tha  conditions  (i)  and  (ii)  and 
let     *|       be  a  giran  rector;  then  there  exists  exactly  one  f\]netion 
with  the  following  three  properties t 


yz 


*•       9U)      is  eontlnuotts  and  continuottsly  differantiabl*  forx*\,a,b], 
b.       7'U)  «|(x,  7^))     ,     )i*La,tl, 

(!•••  th«  initial  Tilu*  problen  giran  b7 

y'=|cx.y)     ,       y^'^-)'  \   >  (76) 

has  a  VBiiqua  solution. )• 

L«t  vs  now  consider  the  RungeJCutta  method  for  oar  systea  of  equatioas. 
¥a  let  X  •  ta.bj    ,  y   be  an  arbitrary  Teetor  and  z(*)    denote  the  solution  of 
the  systea  of  equations  giren  by. 


1'  « 


and  set 

_         ,        I         A 

fu.y)  ^    A  -^0  , 

Wa  can  ^  the  exact  relative  increment  of  the  solution  of     e'  «  4  C-t,^). 
A  one^tep  method  (Range-Kutta  is  the  most  sophisticated  of  these  wthods) 
for  the  solution  of  the  initial  ralue  problem  given  in  equation  (8)  is 
defined  by, 

y-  M    ' 

y-  =     y"    4-    >(»  •$>(<,  y"  i A)    ,     w«o,t,..,    . 

Here  ^  is  called  the  increment  function  and  is  chosen  so  as  to  approximate 
£k     as  dosoly  as  possible. 
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In  ord«r  to  •llalnaU  the  special  role  played  by  the  Independent  vari- 
able x«  «•  augaent  the  system  of  equations  (69)  by  the  differential  equation 

y"(x)  =  1  (77) 

to  be  saUsfled  by  a  new  function     y'CK) ,  subject  to  the  Initial  condition 

y*(x.')  =  X.  (78) 

Equations  (77)  *nd  (78)  clearly  imply  that  y"(x)  «  x.  Hence  we  can  replace 
the  system  of  equations  (69)  by  the  equivalent  system  for  (s-t-1)  fonotians, 

*•"  rcr.y v*)  =  i 

The  system  now  given  in  equation  (79)  has  the  advantage  that  the  vari. 
ables  entering  into  the  functions  \     may  all  be  considered  dependent.  For 
sii^plicity  we  can  continue  to  denote  the  dependent  variables  by  V  ,y\...,y' 
whether  or  not  one  of  them  is  x.  Thus  we  may  write  the  initial  value  pro- 
blem as 

y'^{(y)   ,  Y<^.>^  7  »  (80) 

where  y  ,  t  and  ^  are  all  vectors  with  s  components. 

If  I  (y)  does  not  depend  explicitly  on  x,  neither  does  the  function  ^  , 
aer  does  the  increment  $  .  Hence 

Mew  if  the  f^tuiction  y  (x)  is  a  solution  of  equation  (80)  and  assuming  ■ 
the  oomponents  of  f  are  sufficiently  differentiable  th«i  the  higher  deriva- 


3* 


tlTM  tf  ^(.x)  c«n  be  •xpressed  In  terms  of  the  funotion  -I-  and  its  derlTa* 
tires*  For  exjwpLe 

J  ■  ■«  .  ,        '*.'■".  i"  •. 

Qenerally,  fbr  k  ■  1,2 assuming  differentiability. 


aiBoe  the  Amotions      k    (J)       exist,  we  hare 


(81) 


Aa   in  the  ease  of  a  single  differential  equation,  a  method  for  approxii 
■ate  integration  oan  be  based  on  a  truncated  Taylor  series,  the  inoreaent 
fwotiem  for  the  Taylor  series  expansion  of  order  p  is 

This  method,  being  of  no  great  praotLoal  Interest  since  the  eraluation 
of  many  deilratires  is  InvolTed,  is  hoversr  of  theoretical  importance  slnoe 
the  Rnnge-Xutta  methods  are  based  on  the  idea  of  approximating  equation  (82) 
by  expressions  which  do  not  involTs  any  functions  other  than  \i^)    • 

Iiot  us  look  now  in  more  detail  at  the  deriratlTes  and  their  stnieture 
for   4  (y;  .  To  slaqxLlfy  the  notation,  write  for  i,J,k  ■  1,2,  •••  ,8 
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vh«r«b7   fi  ,  tjk  ,  ...     .  i»  Mant  th©  vectoM  with  the  ooHponents 
f*   r  ,  (1  ■  1,2,  .  .  .  ,8),  Wb  also  will  find  it  appropriat* 

at  this  tlM  to  adopt  tho  suBunatlon  convention  whloh  Is  used  quite  often  in 
Teotor  and  tensor  analysis  (i.e.  if  an  index  occurs  both  ^b   a  subscript  and 
a  supersoript,  the  teras  should  be  summed  with  respect  to  this  Index  from 
1  to  &)•  Thus  we  hare. 

It  is  clear  that  sums  of  products  of  the  font  of  equation  (84)  can  bo 
dirferiQtiated  like  ordinary  products.  Hence 

ut 

where  i.j.k.n  «  1,2,   •••  ,b  and  the  argument  of  eyery  fluiction  is  understood 
t«  by    y    •     Also  denote  by    A  ,B,  ...     the  vectors  with  the  ooi^nnents 
A*",  Bi*",   •••   (1  ■  1,2,   ,,.   ,s)  respectively. 

With  these  oonventions,  we  can  write  formally  the  first  few  derlvatlTSS 
in  the  following  co^iact  manner: 


■V 
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aino*  v»  will  hav»  to  wcpwid  •xpresslons  of  the  font      4  (  Y  ♦  -^a )        In  powers 
of  h,  whoro    y    and     oL    aro  fixed  vectors,  we  write  Taylor's  series  for  func- 
tions of  sereniL  rariables     [16.227]   . 


Vb  want  to  oonbine  values  of  the  ftmction  which    I    takes  at  different  points, 
in  such  a  way  that  the  resulting  function     1^<7,1\)   agrees  as  closely  as        . 
possiULe  with 

(HoU  that  equaUon  (88)  is  the  sam  as  equation  (81).) 
For  the  second  order  Runge-Kutta  aethod,  we  put 


where  a^,  a,  end  p  nust  be  detendned.    Using  equation  (8?),  ve  hare, 

.  =    ;[  *  ^pB  *    y:£)*t    *   O(i')  . 


B«oe 


^(,,1.)  -  u,-..)A  ♦  vi^s  .  -^..ttfV?. »  ou') .-;,. 


Mow  equating  the  constant  and  linear  tent  in  h  (it  is  iaqjosslble  to 
obtain  acreenant  In  h    since     D      is  not  present  in  the  aboTS  equation  )Mith 
tho  oorresponding  tens  in  equation  (88),  we  have. 
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siting  th«  g«n«ral  soltttion  of  tli*  abov*  systea  of  oquAtlons  *s 

vhor*     <A  :^  o    . 

Thus  Dm  iBoroBont  funotion  is  giTon  bj, 

?()>»i^)   *   0-'<)r(y)    *   Mf(y    *    ^  ?Cy))    ,    »<  ^  o    . 


.1 


.  .-'?.•; 


(89) 


It  doriatM  froa  •qtutlon  (88)  by       OdO*)     an'd  two  ovalufttions  of   \{9) 
aro  r«q«lr»d  to  oonputo      ^     • 

la  a  siailar  aumar,   as  In  the  precoedlng  casa  and  analogous  to  tha 
•inglo  aquation,  wa  vill  arrira  at  tha  classical  RungaJCutta  fourth  ordar 
■athod. 

Tha  laoraaant  fnaotlon  Is  glTon  by 

and  ona  aat  of  appropriata  choloas  of  a^ ,  a  ,  a^  and  a.  is  found  to  ba 

'3 


*1  ■%■  1/6     ,     a^  «  a,  ■  1/3     , 


Tins, 


^(r»  =   i(k.^  2^.  *  "2.1;.  I J    ^ 


(90) 
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It  la  d.««r  that  •quatlon  (90)  is  a  special  eas*  of  th«  oLassleal  Rung*- 
IvttA  fourth  oi^or  equation  dorlTod  earlier  for  the  ease  of  the  single 
dlffarmtlal  equation.  >:• 
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NUMERICAL  E3CAMPLES 

Fbr  the  sake  of  Illustration,  we  list  five  n\u»rlcal  exaaqples  (Tattle  1). 
In  papers  of  this  kind,  numerical  examples  are  desirable,   especially  those 
lihleh  Illustrate  how  well  the  derived  method  compares  with  others.    It  is 
sometimes  difficult  to  choose  meaningful  examples  to  Illustrate  Runge-Kutt* 
methods,   since  the  oomplieated  nature  of  the  error  term  makes  it  difficult  to 
dKMse  a  function  f(x,y)  which  really  serves  as  a  test  while  at  the  same  time 
yields  a  problem  which  can  be  solved  analytically. 

A  rOKTRAN  program  was  written  to  con^mte  both  the  classical  and  the  min- 
imum  fourth-order  solution  for  a  given  differential  equation  and  to  calculate 
the  error  in  eadi  method  for  the  desired  number  of  iterations  (see  Appendix). 
The  results  for  these  five  differential  equations  are  given  in  Table  2. 

The  first  three  of  these  exaiqjles  show  that  the  minimum  method  conqiares 
favorably  with  the  classical  method  while  in  the  fourth  exai^ple,  there  is 
really  no  comparison  and  finally  the  fifth  exaB?)le  is  not  nearly  as  favorable 
in  oon)arison.    We  note  that  it  is  only  a  matter  of  a  little  ingenuity  to 
find  other  exaaples  to  make  the  minimum  method  appear  more  or  less  favorable 
in  oo^>arlson  with  the  classical  method  or  other  methods. 

I' 

In  oondosion,  we  re.state  the  main  point  of  this  report.  If  Runge- 
Kutta  methods  are  to  be  used  to  start  the  solution  and/or  to  change  the  inter, 
val  site,  one  is  interested  only  in  being  able  to  limit  the  truncation  error 
to  as  small  a  quantity  as  is  possible.  Hence  we  dioose  the  method  which 
puts  the  smallest  bound  on  the  error  term  in  this  sense.  Therefore,  when  « 
fourth-order  method  is  desired,  equation  (68)  should  be  used,  vhen  a  third-  : 
order  method  is  employed  equation  (61)  should  be  used,  and  equation  (30)   V": 
should  be  used  when  a  seoond-order  method  is  under  consideration.  In  all 


v-M 


itO 


««0M,  if  th«  iMthod  is  best  for  a  single  equation,  it  is  at  least  neaxly 
best  for  a  systen  of  equations. 
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Kxaapl* 


n 


HI 


IV 


TABLE  1 


Differential  Equation 


g  ■  iktlljt-Si 


dx 


d£ 


^    JStj 


dx       1  +  tan*^  y 


dx 


1-/ 


S  -  3y2(xe^  .  6) 


Initial 
Condition 


yd)  «1 
y(0)  =    -1 
y(0)   »  0 
y(0)  3  0 
y(o)  =  1 


Solution 


■  xT«g  X  +  2x      .X 


1  -  W^  »  2x 

7  "  ^k 


y  «  arotan  x 


y  s  tanh  x 


6  -x»* 


1/3 


In  each  of  the  aboTs  problems,  we  calculate  the  value  of  y  vfaen  x  »  k. 
We  calculate  the  "exact"  value  from  the  solution  given  in  Table  1   and  then 
calculate  the  "approximate"  solution  using  numerical  methods.     Thus  In  ex« 
ai^ile  I  ve  calculate  y  for  x  »  1(.1)^  and  also  for  x  «  iCi^^'f  •     In  exampl* 
H  ve  calculate  y  for  x  »  1(.1)4,  etc.    The  results  for  x  *  ^  are  ooq)ared 
vlth  the  "exact*  values  and  these  differences  are  noted  in  Table  2. 


*  Here  the  notation  x  >  iCO'f  means  that  ve  let  x  take  on  successive  values 
tluroughout  the  interval  Arom  XBltox«4ln  steps  of  length  .1 
(i.e.  Xq  ■  1.  x^  ■  1.1,  xg  «  1.2,   ...,  x^^^  «  3.9,  r^  «  4.0). 


kZ 


TABLE  2 


fiUBpl* 

St«p 

Number  of 
Iterations 

Error  for  daasical 
Fourth-Order  Method 

Error  for  Mlninnui 
Fourth-Order  Method 

I 

.1 
.2 

20 

-.338  X  10-3 
-.433  X  10-2 

-.265  X  10-3 
-.328  X  10-2 

n 

.1 
.2 

20 

.200  X  10-^ 
-.900  X  10-^ 

.000  X  10-99 
-.900  X  10-^ 

ni 

.1 
.2 

20 

-.800  X  10-^ 
-.180  X  10-5 

-.700  X  10-^ 
-.150  X  10"5 

IV 

.1 
.2 

i<0 
20 

..000  X  10-99 
.000  X  10-99 

.000  X  10-99 
.000  X  10-99 

V 

.1 

10 

.^^o  X  10-5 

,260  X  10-5 

*3 
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APPENDIX 


KSU    l'»lO    COMPUTING    CLNTLR 


0000  I 

00002 
0000  J 
0000<» 
OOOO*} 
00006 

00007 
OOOOP 
00009 
OOOIU 


DIMtNSItlNlJ(Al) 
FU<MAU  IH     ,7X, 

L.rt) 

FO-tMAT(  IHT,9X, 
rU'<M'\r(  IHS,7X, 
FU-tMAI  (  IH     , 
FD-^MAT  t  IHl, 


I7HARS0LUTB 


47 


ERROR    =    El<f.8,5X,17HRELATIVE    ERROR    =    EH 


bX, 


I 


00018 

000?? 
000?  <► 
000Z6 

0002<» 
00031 

00036 
000*0 

000A5 
00047 


00052 
000^3 


I  MmUJE    SOLUTION) 

RHX    VALULS8X,RHY    VALUES) 

r  lO.B»5X,hlb.8) 
12X,^0HRUNGE-KUTTA    SOLUTION//)  ctcocv     i -.u  f  laTcuMcr 

rcHM\TllHS,7X,(*HX    VALUFSRX.BHY    VALUES5X, 9HSTEP    SIZE25X ,  12H INTERMEI 
lAltlOX.RHX    VALUFS5X,8HY    VALUES) 
FU^MAKIM    ,3X,FI"->.8,5X,F15.8,10X,F3.2) 

RV^K!l[(lHi;4Xjr7HKALiT0Ns'HiN!MUM    RUNGE-KUTTA    SOLUTION//) 

gVv'?xIy)  =  1?/?1.*(SIN(Y)/C0S(Y))»(SIN«Y)/C0SIY))) 

WKIIt(3,?) 
Wl<  I  lk(  3«  3) 

0i)U'K=l,40  ■" ~ 

AN-K 

U(K)=AN«.l 

HtCv  )^ATAN(U(K)  ) 

WKlTFil  3,<»)U(H)  ,W(K)  ' 

CONTINUE  -^c' 

W><  till  3,6) 

Ri^  A.J(  l,l0)XA,YA,DA,XAND,YAND 

IF  (4)A.EO.O.)GOT060 

DO     .H»DA 

X*  <A 
Y  =  YA 
J=  I 
/=  .IV( 


■II 


U' 


X,Y) 


00060 


OOOfeJ 
000?.4 

00066 


00069 
00071 


r,ulu(3lf  36f40t<.^),J 

CI -u^•l 
x=xA»on 

Y=YA»Cl/2. 

G0IU29 
C2-UA«Z 
Y=YA*C2/2. 
J  =  3 

GOTU?g 
C3=UA»Z 
X  =  XA«-OA 
Y=YA*C3 

G0r029 

XA  =  X 

C'^  =  DA«Z 

YA  =  Y^MCl♦2.•(C2♦C3)♦C4)/6. 

JJ  -JJ*1 

IFt JJ.LT.25)G0T052 

WRI  TEI 3,8)XA,YA 

JJ  =  0 

IF(XA.LT.XAND)GUTU26 

XANO=XA 

YAND=YA 

ABrK=YAND-W(40) 

R6Lb=AUER/W(40) 

WRI TEI 3,7)XANU,YAND,DA 

WRITE(3t  DABEKfRELE 

G0T(J?2 

CllNTlNUr 

WK  Mt  ( 3t9) 

W«  MCI  3,6) 

Rt  Aul  l,l0)XB,YB,Dlif  XBNDfYBND 

IF(UB.EQ.O.)STUP 

KK  =  0 

X  =  Xb 

Y  =  YD 

I'l 

Z=GIVIX,Y) 

G0ru(71(76»81t86),I 

Cl=OB»Z 


KSU    I^IO    Cl)MPUTING    CTNTER 

Guruis^  -. 

00076    C«?  =  Ut^«Z 

X  =  XU*.'.5^7:i7»OB 

Y  =  YH*.2969  78»CU.lb8759»C2 

1=  ) 

GUTUb9 
00081    C3--u;i»Z 

w  _  y  11  A.  i  )  Ll 

Y=YtU.2l8l00»Cl-3.050965»C2+3.832864»C3 

Guru69 

00OP6    Xl^^X 

00088    YR=Ya"f.l7^760»Cl-.55l481»C2+1.205536»C3*.l71185»C*) 

IFUK.lt. 25)G0TU93 

WrtlTFC 3,8)XB,YB 

KK  =  0 
00093    irt Xfl.LT.XRMD) G0TO66 
0009*    XONU^XB 

YBNl)=YR 

APbK=YBNO-W(*0) 

RELE  =  ABf:R/W(*0) 

WHlITEl  3,7)XBND,YBND,DD 

WRI r6(3tl)ABER,KELE 

G0TU63 

END 
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Runge-KutU  Methods  are  numerical  means  for  solving  a  differential 
equation  (or  a  system  of  differential  equations)  given  Initial  values.     They 
are  one^tep  methods  and  require   (for  the  most  popular  method)  four  Itera* 
tlons  at  each  stage  of  the  calculation.     Con?)ared  with  the  more  popular 
■ulti-step  methods  which  require  two  iterations  for  each  calculation,   Runge- 
Kutta  methods  are  more  time  consuming  even  on  present  day  coBqmters.     Runge- 
Kutta  methods  are,  however,   self-starting  and  as  andx  are  \ised  primarily  to 
calculate  starting  values  to  be  used  then  by  the  more  stable  predictor, 
corrector  or  multl.step  methods. 

Considering,  then,   Runge-Kutta  methods  only  for  starting  the  solution, 
we  are  concerned  with  being  able  to  mlnlmlae  the  truncation  or  discretiza- 
tion error.     In  this  report,  we  derive  Runge-Kutta  methods  of  second,  third 
and  fourth  orders,   and  use  this  derivation,   assuming  certain  bounds  on  the 
function  and  its  partial  derivatives,  to  arrive  at  expressions  for  the  error 
term  in  each  method.     These  are  minimized  by  appropriate  choices  of  the 
arbitrary  parameters.     With  these  choices  for  the  arbitrary  parameters,  new 
ooeffioiants  are  determined  and  used  to  write  the  minimum  Runge-Kutta  methods. 

in  analogous  treatment  of  the  derivation  for  a  single  differential 
•qoation  is  given  for  a  system  of  differential  equations. 

Five  differential  equations  are  chosen  as  examples.     Each  differential 
equation  was  solved  by  the  classical  method  and  by  the  minimum  method,   and 
the  erzvr  was  calculated  in  each  case  after  a  particular  number  of  itera- 
tions.    A  FORTRAN  program  was  written  and  the  I'flO  computer  was  used  to  carzy 
out  the  computations.     Although  some  of  the  examples  conqiare  more  favorably 
with  the  theory  than  others,  it  was  pointed  out  that  with  a  little  foresight, 
one  oan  find  additional  e(xai^>les  whidi  either  do  or  do  not  ootapT9  favorably 
with  tha  UMory. 


If  Runge-Kutt*  methods  ar«  to  be  used  to  start  a  solution,   these  Mini  mm 
methods  are  generally  better.     Jttso  if  a  system  of  differential  eqnaUona  is 
under  oonside ration,  these  methods  will  be  at  least  nearly  best  In  the  mini, 
mum  truncation  error  sense. 


